Numerical simulation of oscillatons: extracting the radiating tail 
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Spherically symmetric, time-periodic oscillatons ~ solutions of the Einstein-Klein-Gordon system 
(^— N \ (a massive scalar field coupled to gravity) with a spatially localized core - are investigated by 

fvj ■ very precise numerical techniques based on spectral methods. In particular the amplitude of their 

' standing-wave tail is determined. It is found that the amplitude of the oscillating tail is very 

OJJ' small, but non-vanishing for the range of frequencies considered. It follows that exactly time- 

periodic oscillatons are not truly localized, and they can be pictured loosely as consisting of a well 
(exponentially) localized nonsingular core and an oscillating tail making the total mass infinite. 
Finite mass physical oscillatons with a well localized core - solutions of the Cauchy-problem with 
suitable initial conditions - are only approximately time-periodic. They are continuously losing their 
mass because the scalar field radiates to infinity. Their core and radiative tail is well approximated 
by that of time-periodic oscillatons. Moreover the mass loss rate of physical oscillatons is estimated 
O , from the numerical data and a semi-empirical formula is deduced. The numerical results are in 

<vJ i agreement with those obtained analytically in the limit of small amplitude time-periodic oscillatons. 
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PACS numbers: 02.70.Hm, 03.50.Kk, 04.25.dc, 04.40.Nr 



I. INTRODUCTION 



0^ ■ Oscillating soliton stars with spherical symmetry were first observed in numerical simulations performed by Seidel 

C ' and Suen in [1]. They considered a free, real massive scalar field (of mass m) coupled to gravity. Starting from general 
^^ initial data, it was found that the system could settle in a spatially localized and apparently time-periodic and stable 
t^^ ' state. Such spatially localized, oscillating solutions are referred to as oscillatons. There seems to be at least a one 
^D [ parameter family of oscillatons, labeled by their frequency, lu < m. The numerical simulations were consistent with 
^~~^ • oscillatons being time-periodic, with a constant frequency. Assuming exact periodicity the structure of such solutions 
\ I ' was then investigated by means of a Fourier decomposition of the various fields describing the system. At that time, 
^ it was not clear if oscillatons were truly periodic and Seidel and Suen left open the possibility of the appearance of 
a secular change of their frequency. The same authors studied the formation of oscillatons in a scenario involving 
^ . , gravitational cooling in the companion paper [2], making them of possibly great physical importance. Oscillatons 
j^ ' appear to be good candidates for dark matter in our Universe [3-8] . 

Some years later, this subject was re- investigated in a series of papers [9-11]. Several aspects of oscillatons were 
studied, like their structure, dynamics, stability and more. In particular it was found that oscillatons reached a 
maximum mass around an oscillation frequency of Wmin ~ 0.86m. It is expected that solutions would become unstable 
below this frequency, therefore the frequency range of oscillatons is expected to be Wmin < uj < m. The possible 
existence of a secular change of the frequency of oscillatons induced by (scalar) radiation has not been discussed, 
however, in those papers, as a matter of fact it has been simply assumed that the solutions are truly periodic. 

The physically motivated, important question concerning the possible radiation loss of oscillatons, has been first 
considered by Don N. Page [12], who has pointed out that due to scalar radiation oscillatons are bound to lose mass. 
Page has estimated the amplitude of the outgoing scalar wave and derived a formula for their mass-loss, which turns 
out to be rather small even on cosmological time-scales. The subject of longevity of oscillatons has been taken up 
again in a recent paper [13]. This study has also been performed in the limit where the amplitude of the scalar field is 
small, when a perturbative analytic approach is feasible, and the structure of oscillatons was investigated in detail. In 
the limit of small amplitudes, oscillatons correspond to solutions with frequencies a; — >■ m. The fact that oscillatons 
are continuously losing mass due to scalar radiation, leading to secular changes in their frequency has been confirmed. 
A formula for the mass loss of oscillatons in spatial dimensions 2 < D < 6 has been given. Although the quantitative 
result of Ref. [13] for the mass loss rate differs significantly from that of Ref. [12], there is qualitative agreement on 



its extreme smallness as compared to the total energy of the configurations. This explains why physical oscillatons, 
evolving from finite energy initial data, are practically indistinguishable from exactly time-periodic oscillatons, which 
have an extremely small amplitude standing wave tail, and also why this effect has been largely ignored by previous 
numerical work. 

In various field theories containing massive scalars, very similar objects, oscillons or pulsons have been numerically 
observed and analytically investigated mostly in the limit of small amplitudes [14-25]. Oscillatons under the name 
of "gravipulsons" have been studied in the recent paper [26] in EKG theory with a logarithmic scalar self-interaction 
analytically without assuming a small amplitude limit 

The goal of this paper is to investigate exactly time-periodic oscillatons in 3 dimensional EKG theory and generalize 
some of the results of Ref. [13] without using the small amplitude limit. The asymptotic oscillating tails of the fields 
are determined by special high precision numerical techniques. The extreme smallness of the oscillatory tail of an 
oscillaton with respect to its central amplitude [13] explains why it has been missed in previous numerical work. 
We show that nevertheless this tail can be computed by making use of very precise numerical techniques. In this 
context, the field equations are solved using spectral methods, where the variables are approximated as finite sums 
of known functions called the basis functions. This is done for both space and time. The solution of the numerical 
system is performed by using the spectral solver Kadath [27, 28] which enables the use of spectral methods in a 
wide range of problems arising in theoretical physics. The oscillating tails of the scalar fields are obtained by using 
a similar approach as the one used in Ref. [29] in the case of oscillons in a self-interacting single scalar field theory 
in 3 dimensions. A detailed comparison of the numerical results with those obtained in the small amplitude limit 
shows remarkable agreement and coherence. The core of the oscillaton is approximated by the expansion in the 
small amplitude limit to good precision even for not too small values of the amplitude. This expansion significantly 
underestimates, however, the magnitude of the oscillating tail which depends in an essentially non- analytic way from 
the amplitude. 

The paper is organized as follows. In Sec. II the equations describing a scalar field coupled to gravity are pre- 
sented, along with the decomposition in modes used to obtain periodic and weakly localized solutions. The expected 
asymptotic behaviors of the various fields are also studied. Section III is devoted to the presentation of the numerical 
techniques. The use of the library Kadath and the way solutions are matched at the outer end of the computational 
domain are explained in some detail. Numerical results are shown in Sec. IV along with many numerical tests that 
validate the overall procedure. The existence of an oscillatory tail and comparison with previous work are discussed. 
The comparison with results from the small-amplitude expansion is shown in Sec. V. The mass loss rate of oscillatons 
is determined in Sec. VI. 

II. THE MODEL 
A. The equations 

One considers a real scalar field $, coupled to gravity. The stress-energy tensor is given by 



r^^ = ^.M^.-^ - 9i^iy 



i$ „$'" + [/($) 



(1) 



where C/($) is the potential of self-interaction. 

We are interested in finding spherically symmetric configurations in 3+1 dimensions. We can chose the metric to 
take the following form : 

ds^ = -Adt^ + B (dr^ + r^dn) , (2) 

where fl is the solid angle. The unknown functions A and B depend solely on r and t, and this system of coordinates 
may be described as quasi-isotropic. 

In order to get rid of the Stt factors (see Sec. II- A of [13]) we rescale the scalar field and the potential as 

<i>^^a>, u{<P)^^u{^). (3) 
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Under these assumptions, Einstein's equations G^,y = SttT^^ are written as : 
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Equations (4, 5, 6, 7) correspond respectively to the components 2Ett, 2Err, Etr and 2 [Eee/r^ — Err) of Einstein's 
equations. Equation (8) is the conservation of the energy-momentum tensor. 

In the foUowing, one wih consider the simplest choice for the potential, U — ?ti<I>2/2, a free massive Klein-Gordon 
field, and by suitable scaling we set to = 1. 

B. Mode decomposition 

As in Ref. [1], one wishes to study the possibility that there exist periodic, or quasi-periodic solutions, to the system 
(4-8). Following [1, 13] one assumes that A and B contain only even harmonics with respect to time and $ only odd 
ones, which is valid for any symmetric potential U . This gives the following mode-decompositions : 



A(r,t) 

B{r,t) 

$ (r, t) 
where w is the frequency of the solution. 



l + 5I^2j Wcos((2j>t) 



3=0 



l + II^2j(r)cos((2i)wt) 



3=0 



J2 *2j+i (r) cos ((2j + 1) ujt) , 



(9) 
(10) 
(11) 



C. Asymptotic behavior 



The asymptotic behavior of time periodic solutions of the EKG system is complicated by the fact that the existence 
of standing wave tails is not compatible with asymptotic flatness. Heuristically, no matter how small the amplitude of 
the oscillating tail, due to its slow spatial decay the total mass is infinite. Fortunately one can sidestep this somewhat 
complicated issue, since there is an intermediate asymptotic region, which is defined by being sufficiently far from the 
core region of the oscillaton where the first Fourier component of the oscillating tail still dominates, but the mass in 
this tail is negligible with respect to the mass in the core. In Ref. [29], we succeeded to determine numerically the 
amplitude of the standing wave tail of oscillons in the case of a single scalar field with a self-interaction potential in 
fiat space-time. In this context the amplitude of the oscillons was successfully isolated using a formalism based in 
the homogeneous solutions of the various operators. The aim of this work is to use similar techniques in the case of 
the oscillaton. The first step is to assume that there is an intermediate region where space-time can be considered 
asymptotically flat. This translates to the following behavior for the metric fields, at large radius : 

(12) 

(13) 

For a solution of the full system, one must have ta = tb = fQ. This is however not enforced directly in the numerical 
solution but rather used as a measure of the accuracy of the code. 
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The situation of the scalar field is more complicated and one has to study in some detail the wave-equation (8). If 
one keeps only the dominating terms (i.e. if one sets A = 1 and B = 1), Eq. (8) reduces to : 

$,rr+-$,r-$.U-* = 0. (14) 

r 
$ being a sum of odd cosines, the equation for the harmonic $„ is : 

(r$„)_„ + {nW - 1) r$„ = . (15) 

Given that w < 1, one has two different cases : 

• For n = 1, the solution that vanishes at infinity is 

$r(.) = C,2^^Pi^, (16) 



with e = Vl — t^^- 
• For n > 1, the solution is oscillatory of the form 



-i>„(r)^c /°^^'"^ + ""\ (17) 



with A„ = \/n?uj^ — 1. 



At this order of approximation, the background is fiat Minkowskian, and the phase a„ is constant. However, 
considering waves on a Schwarzschild background, the phase will have a slow radial dependence. Hence, it will prove 
useful to get the next order of approximation for n > 1. In order to do so, we consider the phases a„ as slowly varying 
functions of the radius r, so that an,r ^ A„. More precisely, one assumes that an,r is of order A„/r. Keeping the 
first two orders in terms of 1/r, one gets : 

^ / \ ^ ^^ , , sin(A„r + a„) cos(A„r + Q;„) 

$«,r(r) == -C„(A„ + a„,r) ^ --Cn ^ (18) 

rt, t \ n (\^ ^o\ \ cos(A^^r + a^ sin (A„r + a„) 

So it appears that the terms involving «„ ,. are of the same order as the corrections induced by the r^/r parts of the 
metric fields, the equation for the scalar field, up to 1/r^ terms, being 

1 - -) $«,rr + -*„,r + n^W^ (^^ + I^^ $^ _ $^^ = Q . (20) 

When inserting Eqs. (18-19) into (20), the second order terms lead to an equation for the phase «„ (the first order 
condition has already been used and gives the value of A„). One finds that 

-2A„a„,, + ^ (A2 + n^cj^) = 0, (21) 



which can be integrated to give : 



logr + ,5„, (22) 



where 5n is a constant. One can check that, as assumed, a„.r is indeed of the order of A„/r, in agreement with the 
starting hypothesis. A detailed study of series solutions for the Klein-Gordon equation on Schwarzschild spacetime 
can be found in [31]. 

The oscillatory behavior of Eq. (17) illustrates the fact that, in general, time periodic solutions cannot be truly 
localized, nor can the corresponding space-time be asymptotically fiat. Indeed, if $„ behaves like (17), the terms 
involving $ in Eqs. (4-7) will decrease only like l/r^, thus being inconsistent with the dominating behaviors of ^ and 
B assumed from the start and given by Eqs. (12-13). It is only due to the smallness of the amplitude of the tail that 
an "intermediate asymptotic" region exists. 

In some very particular situations, it can occur that the oscillatory tail is absent thus truly localized time periodic 
solutions do exist. A famous example is the breather solution in sine-Gordon theory in 1 dimension. Results from 
[1, 2, 13] strongly suggest that, even if the oscillatory tail is present, its amplitude should be relatively small, so that 
there exists a region where its influence on the various fields is also small. 



III. NUMERICAL METHODS 
A. Spectral expansion 

Solutions of the system are sought by making use of the spectral library Kadath [27, 28]. The setting uses a 
two-dimensional space, with respect to the coordinates {t,r). For the time-coordinate, a single domain is used. The 
physical time t relates to the numerical one t* hy t* ^ uit. A spectral expansion is then performed with respect to t*, 
using only even cosines for A and B and only odd ones for $, in accordance with the mode decompositions (9-11). 

For the radial coordinate a multi-domain decomposition is used, similar to the one described in Sec. (2.2) of 
[28] . Typically one considers a nucleus that contains the origin and several spherical shells that are bounded by two 
finite radii. However, given the appearance of oscillatory solutions (see Sec. 11 C), no compactification of space is 
used and the equations are solved only up to a given radius i?max at which an appropriate matching is performed 
(see Sec. TUB). In each domain the physical radius r is related to the numerical one r* by an affine-law. In the 
nucleus one uses r = i?nuc x '''* with r* S [0, 1] and where i?nuc is the radius of the nucleus. In the shells one uses 

-'^outcr -'^^inncr \ ^ / -'^outcr ~'~ -'•'inner \ -iiTt^r-i-ii 11 t~, in ,1- i , 

With r G [^1,1] and where rtinnor and rtoutcr are the inner and outer 
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radii of the domain. Spectral expansion is performed with respect to r* . In the nucleus, and to account for the fact 
that the fields are even near the origin, only even Chebyshev polynomials are used (this is also the reason for the 
different range of variation of r*). Spectral expansion is performed with respect to standard Chebyshev polynomials 
in the various shells. 

For instance, in a given shell, A is approximated by A « >, /, ^ij cos {2jt*) Ti (r*). Nt and A^^ are the number of 

i=i i=o 
coefficients with respect to t* and r*. Ti denotes the i*"^ Chebyshev polynomial. The Aij are the spectral coefficients 
of ^. 

Let us finally point out that all the divisions by r that appear in the equations concern quantities that are odd near 
the origin (like A^r)- Such ratios are then easily computed using the spectral expansion because the ratio of T2i+i/r 
can be exactly expressed as a sum of T2i. The division is said to take place in the coeflficient space. 

B. Matching criteria at the outer boundary 

Given that the computational domain can not be extended to infinity, one needs to derive appropriate outer 
boundary conditions for the various fields. They are based on the asymptotic behaviors in the "intermediate region" 
found in Sec. II C. The idea is to match the various fields to the solutions of the dominating operators, at large 
radius. More explicitly, let us consider a purely radial function / (r) that must be matched to the solution g (r) of the 
operator. By hypothesis, one demands that, at a large matching radius R, f is close to a solution of the form Cg (r), 
where C is a constant not known a priori. The continuity of / and its radial derivative f^r at the matching radius 
gives the following conditions: 

/ (-Rmax) — Cg (i?max) (23) 

f,r (-Rmax) = Cg^r (-Rmax) • (24) 

The constant C can be eliminated to get a boundary condition for / that involves only the function g : 

[fg.r - f.rg] (R) = 0. (25) 

Once the function / is known, the constant C can be recovered by making use of either Eq. (23) or (24). 

This matching technique is used to match A (resp. B) to a solutions of the type 1 (resp. 1 ), where 

r r 

the constant rA (resp. r^) plays the role of C in Eqs. (23) and (24). Let us note that for those two fields, one only 

matches the harmonics which are time-independent, the other ones being simply set to zero. 

For the scalar field, the matching functions depend on the order n of the harmonic considered and are given by Eq. 

(16) for n = 1 and Eq. (17) otherwise. In the expression (22) of the slowly varying phase, one will use rs for the 

value ro, as both rA and rs are expected to converge to rg as the matching radius i?max increases. 



C. Minimization of the oscillatory tail 

In the asymptotics given in Sec. II C, for each n > 1, there is a constant phase 5„, that can be chosen freely (see 
Eq. (22)). In practice, one wishes to construct solutions that are as close as possible to truly localized ones. Therefore 
it is desirable to make the oscillatory tails as small as possible. As seen in Sec. II C, the dominating mode for which 
the oscillations appear is $3 therefore we try to minimize C3 (C, being the amplitude of the oscillatory tail of the i*"^ 
mode, as seen in Sec. IIIB). 

An additional simplification comes from the fact that the oscillatory tails of the other modes n > 3 have a very 
small influence on C3, since their amplitudes, C„, are decreasing very rapidly for increasing values of n. In other 
words, C3 is almost independent of (5„ for n > 3. So the minimization of C3 as a function of only one constant phase 
63 is a very good approximation to the true minimum. This is done using a standard golden section search algorithm 
(see for instance [30]). The amplitude, C3, never gets so small that the other C„'s need be considered (see Sec. IV C). 

D. The numerical system 

In the context of the spectral methods implemented by the library Kadath, a system of partial differential equations 
on the fields is transformed into a set of algebraic equations for the unknowns that are the coefficients of the various 
fields. The non-linear system is solved by making use of a Newton-Raphson iteration: starting from an initial guess 
that is as good as possible, the solution is found by iteration. At each step, the linearized system (with respect to the 
unknowns) is inverted (see Sec. 5 of [28] for more details). 

The set of equations (4-8) is redundant, meaning there are more equations than unknowns. In order to select an 
appropriate subset of equations let us recall that the actual number of unknowns of a field depends on the exact 
spectral basis, some coefficients not being true degrees of freedom (see Sec. 3.5 of [28] for a detailed discussion on 
that) . The same is true for the equations that each give a number of algebraic equations that depend on the spectral 
basis of the result. In order to ensure that the number of unknowns is identical to the number of equations, one then 
needs to select equations that have the same spectral basis as the fields. Eq (8) has the same basis as the scalar field 
itself, whereas Eq (6) has both a different radial and temporal basis from the ones of A and B (because the equations 
involve one derivative in each dimension). However the three equations (4,5,7) are consistent with the basis of the 
metric fields. There is no argument to discard one more than the other and one just has to check, after solution that 
the "forgotten" equation is indeed verified. It is found that the set of Eqs. (4) (5) and (8) lead to a numerical solution 
that fulfills the full set of equations (sec Sec. IV B). 

In the context of Kadath, partial differential equations are dealt with by means of a tau-mcthod (see for instance 
[32, 33]). In each domain, the equations are solved by demanding that the coefficients of the residual vanish: in 
a sense, one solves the equations in the coefficient space. Depending on the order of the equations, the conditions 
corresponding to the last coefficients must be relaxed to enforce continuity of the solution and appropriate boundary 
conditions. The number of conditions that must be relaxed is known as the order of the method and is closely related 
to the number of homogeneous solutions of the operators. The functions being periodic, there is no need for any 
boundary conditions with respect to time. 

For the radial coordinate, the situation depends on both the equation and the type of domain. The radial derivative 
of the highest order in Eq. (4) is Brr so that one will consider Eq (4) to be associated with the field B. Brr has two 
homogeneous solutions: C and r. However, in the nucleus, the function r is not admissible because, by construction, 
we restricted ourselves to functions B that are symmetric near the origin. So we have only one admissible homogeneous 
solution in the nucleus and two in the various shells so that we have to use a first order tau-mcthod in the nucleus 
and a second order one in the different shells. Both the variable B and its radial derivative Bj. must be matched 
at the interface in the various domains and one must also supply an outer boundary. One can easily check that the 
number of conditions discarded by the tau method is equal to the number of matching and boundary conditions. 

Given that Eq. (4) is associated with the variable B, one has to consider Eq. (5) as an equation for A. The 
highest order derivative is A^r which admits only one homogeneous solution. Equation (5) is then solved by using 
a first order tau-method in every domain, supplemented by the matching of A at each interface and an appropriate 
outer boundary condition. The situation for $ is similar to the case of B. Let us point out that associating variables 
with given equations may seem a bit dubious, given that the whole system is coupled, but it does provide very useful 
guideline in constructing a numerical system that is well-posed. The whole situation is summarized in Table I. 



TABLE I: Construction of the numerical system to be inverted. The first two columns show the formal association of each 
unknown to an equation. The order of the iau-methods in both the nucleus and the shells are then given, as well as the 
quantities that must be matched at the boundaries between the domains. The last column gives the function to which each 
harmonic is matched at the outer radius i?max- 
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IV. NUMERICAL RESULTS 



A. Setting 



As observed in [13], the solutions are more and more spatially extended as w gets closer to 1, so that it is convenient 
to work with the variable p = re, where e — \[\^^u? . In terms of p the geometry of the solutions for various uj is 
rather similar. Using this fact as a guideline, we chose the radial domains to be, in terms of r, [0, 1/e], [1/e, 2/e], 
[2/e, 4/e] for the three first ones and [4i/£, 4 (i + 1) /e] afterwards. The size of the last domains remains fixed to 4/e 
because one wishes to be able to resolve the oscillatory tails that appear in $ and so one must ensure that the size 
of the domains does not get too big, with respect to the wavelength of the oscillations. This setting, along with the 
number of radial coefficients used, seems to be satisfactory, as can be seen in Sec. IV C. The outer matching radius 
can be varied by changing the number of domains N^- As a standard value one uses Nd = 20, which corresponds to 
a matching radius of i?max ~ 68/e. 

As far as the number of coefficients are concerned, one uses three different resolutions, a low one with A'^^ — 13 
radial coefficients and Nt = 5 coefficients in time, a medium one with A^^ = 17 and Nt = 7 and a high resolution with 
Nr = 33 and Nt = 9. The system is iterated until the residuals reach a threshold of 10^^ for the low and medium 
resolutions and of 10~^° for the high resolution. The setting of the threshold is somewhat empirical and is related 
to the precision one can expect given a number of coefficients (which is very difficult to assess a priori in the general 
case). A threshold too big will obviously result in a poor precision, no matter what the number of coefficients is. On 
the other hand, if the threshold is too small, it can prevent the code from converging, if the number of coefficients 
is not sufficient. Depending on the cases, the overall precision of the resolution is limited either by the number of 
coefficients or by the value of the threshold (see the discussion about Fig. 1). 

The constant phases Sn are fixed to zero for n ^ 3, whereas ^3 is obtained by minimizing the coefficient C3. The 
search for the minimum is stopped when a precision of 10~^ is achieved on the value of ^3. It will indeed appear that 
even if C3 varies greatly as a function of 63 , the curve is very fiat near the minimum value so that a greater precision 
on ^3 is not required to get a precise value of C3. For instance, for uj = 0.86, this is sufficient to get a relative precision 
on C3 of about 10"'"', which is smaller than other sources of errors. 

As an initial guess for starting the Newton-Raphson iteration, one uses the first order expansion in terms of e given 
in Sec. III-B of [13], for e « 0.1. Once some solutions are known for w close to 1, they can be used as an initial guess 
for computing new configurations, by slowly varying cv. 



Tests 



In order to assess the precision reached by the numerical code, one can check the residual error on the equations 
that have not been used explicitly in finding the solution, i.e. Eqs. (6) and (7). The maximal error measured by the 
highest coefficient of the residuals in shown in Fig. 1, for the three different resolutions and various values of w. In 
general the lower w, the bigger the various modes are so that one needs more coefficients to achieve a given precision 
(this is similar to the case of the Kerr black hole ; see Fig. 9 of [28]). So, for a fixed number of coefficients the expected 
errors coming from the spectral representation should be higher for lower values of w. This is what is observed for the 
lowest resolution until an error of around 10^^ is achieved. This saturation is coming from the value of the threshold 
set to stop the Newton-Raphson iteration. The case of the highest resolution is also easy to understand. The fact 
that the curve is almost flat at a level of 10"^^ illustrates the fact that, in this case, the factor limiting the precision is 



not the number of coefRcients but solely the value of the threshold, set in this case to 10"^*'. The case of the medium 
resolution can seem a bit puzzling. As for the low resolution, the error decreases when uj increases, showing that one is 
dominated by the errors coming from the spectral approximation. However, contrary to what could be expected, the 
error does not show any sign of saturation at the level of the threshold of the Newton-Raphson iteration. The error 
is even getting smaller than in the high resolution case, for high values of w. The reason for that is to be found in 
the precise way the Newton-Raphson iteration proceeds. The iteration is stopped as soon as the error gets below the 
threshold. In general when this happens, the error is just below the threshold, and so of the same order of magnitude 
(this is the case for the highest resolution for instance). However, in the medium resolution, the first value below the 
threshold is actually much lower that the threshold itself. The Newton-Raphson iteration would have stopped at the 
same iteration even for a much smaller threshold. This is why the curve for the medium resolution does not show 
signs of saturation. This is nothing profound and is only coincidental for this particular problem. 

We could have tried to reduce the residual error for the highest resolution by lowering the threshold but this could 
lead to convergence issues coming from round-off errors, the code using only double precision. Moreover, the precision 
reached is sufficient for our purpose, which is the precise extraction of the oscillatory tails (see Sec. IV C). It is 
worth mentioning that the errors shown in Fig. 1 do not measure the accuracy at which the oscillatory tails can be 
extracted. It will indeed appear clearly in Sec. IV C that the results from the high resolution cases are more accurate 
than those from the medium one, especially when the phase ^3 is concerned. Let us finally mention that using Eq. 
(7) in place of Eq. (5) leads to the appearance of spurious solutions that do not fulfill the full set of equations. 




FIG. 1: Maximum error on the full set of equations, as a function of u, for the three different resolutions. 



Figure 2 shows how the numerical results depend on the matching radius -Rmax- It is a very important test that 
validates the matching techniques explained in Sees. II C and IIIB. In the first panel, the masses are shown as a 
function of i?max- In three spatial dimensions, the masses can be read from either A or B and relate to the constants 
ta and r^ by Ma — ta/'^ and Mb = TbI"^ (see Sec. II C for the definition of ta and r^)- As expected. Ma and Mb 
converge to the same value as -Rmax increases. Ma being an increasing function of i?max whereas Mb is a decreasing 
one. The convergence is clearly seen in the second panel, where the relative difference between Ma and Mb is plotted, 
for Lo = 0.86. The curve goes to zero almost exactly like l/r, as expected. The third panel shows the relative error in 
the value of C3 as a function of i?max for three different values of w. More precisely, we show the relative difference 
with the value obtained for our largest Rmax which is 68/e. For the three curves the minimization with respect to Sg, 
is not performed in each case but S3, is fixed to the value found for i?max = 68/e. The matching procedure is validated 
as the relative difference decreases fast as i?max increases. There is a saturation of the relative error at a level of about 
10"''. The fact that the saturation depends on the resolution (as shown by the dashed curved in the case oi uj — 0.86) 
implies that this is an effect of the overall precision of the code ; one is indeed talking of a relative difference of 10""* 
on a quantity which is already very small (10^^ for lj ~ 0.83 and 10~^° for w — 0.92). Those curves indicate that one 
can expect to reach a relative precision of around 10^ on the value of C3. 

Figure 3 shows the influence of the phases (5„ on the value of C3, for lo — 0.86. The flrst panel shows C3 as a 
function of S3 and one can see that the value varies by two orders of magnitude. The value of the minimum found 
by the golden section search algorithm is indicated by the circle. The second panel shows the influence of ^5 on C3. 
More precisely one shows the relative difference IC3 {S5) — C3 (^5 = 0)| /C3 {S5 —0). S3 is fixed to the value found by 
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FIG. 2: The first panel shows the behavior of the masses as a function of i?max- The increasing (resp. decreasing) curves 
correspond to the masses extracted from A (resp. B). The circles correspond to uj = 0.83, the squares to a; = 0.86 and the 
triangles to uj = 0.92. The second panel shows the relative diflference between Ma and Mb, as a function of the radius, for 
UJ = 0.86. The curve decreases as 1/r as can be seen by comparing it with the bold curve, which is exactly 1/r. The third 
panel shows the convergence of C3 as a function of i?max for three different values of u. In the case of a; = 0.86 we also show 
the results for the medium resolution (squares and dashed curve) . 
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minimization with S^ = 0. The other phases are fixed to zero. As expected it is very small, the variations on the 
value of C3 being of the order of 10^^. This validates the fact that one seeks the minimum of C3 by varying only the 
phase S3. 



co=0.86 
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FIG. 3: On the first panel, for cj = 0.86, the value of C3 is given as a function of the phase 83. The circle corresponds to the 
position of the minimum found with the golden section search algorithm. The second panel shows the influence on the next 
constant phase 65, where ^3 is fixed to its value for ^5 — 0. 



Results 



The main result of our simulations is to show that there is no value of uj for which the oscillatory tail vanishes. 
This can be seen on the first panel of Fig. 4, where the minimum value of C3 is plotted, as a function of cj, for the 
three different resolutions. The fact that C3 does not vanish implies that there exist no truly periodic solutions of 
the system. Nevertheless, C3 is very small and thus one can have very long-lived solutions as observed in [1, 2] for 
instance. The constant phase 63 for which the minimum value of C3 is attained is shown in the second panel of Fig. 4. 
It appears that the low resolution is not precise enough to give a good location of the minimum and that the medium 
one seems to be accurate only for low values of lu for which the various modes $„ are bigger. The results from the 
highest resolution solution are very smooth and one can expect them to be accurate on the whole range of w. Let us 
mention that even if the phase is not always very precisely determined, the value of C3 is obtained with a relatively 
good accuracy, the curve being rather flat near the minimum (see the flrst panel of Fig. 3). 

As an illustration, in Fig. 5, we show the mode $3 as a function of r, for to = 0.86. This mode is the first that is 
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FIG. 4: The first panel sho-ws the minimum value Ca, as a function of u, for the three different resolutions. For the same 
configurations, the second panel shows the phase S3 for -which the minimum is attained. 



matched to an oscillatory solution. The different panels sho"w the field in different regions: i) near the origin in the 
first panel, ii) in transition region "where the oscillations begin to dominate in the second panel and iii) the region 
close to the matching radius i?max in the third panel. In this latter case, the circle denotes the value of i?max- For 
r > i?max there is no numerical solution so we plot the analytical formula (17) instead. The numerical solution and 
its analytical continuation arc indistinguishable by eye. 

Some global quantities are sho-wn in the various panels of Fig. 6. The first panel sho"ws the total mass of the system. 
Given the behavior of the masses shown in the first panel of Fig. 2, we define the mass as the mean of the value 
given for A and B, that is {rA + tb) /4. A maximum mass of Mmax = 0.60535 is attained for Wmin = 0.8608. This is 
consistent with the values given in [10, 11] where the authors found Wmin = 0.864 and Mmax = 0.607. The presence 
of this maximum is important because oscillatons with uj < Wmin are unstable. The second panel of Fig. 6 gives the 
value of $1 at the origin. Finally, in the third panel, we show the transition radius i?trans defined at the radius at 
which the oscillatory tail begins to dominate $3. More precisely, we define i?trans as the first radius for which $3 
vanishes. As expected i?trans increases with w. The wiggling of the curve is just an effect of what the phase of the 
oscillatory tail is when it starts to dominate. 

The various modes are shown as a function of the radius in Fig. 7 for the configuration of maximum mass (i.e. for 
UJ = Wjiiin)- The dominating oscillatory term in $3 is clearly visible. The oscillations that can be seen on the metric 
fields come solely from the coupling with the scalar field. Indeed, remember that no matching with oscillatory tails 
is used for the metric fields, as explained in Sec. IIIB. When we tried comparing the values of <&„ with the plots of 
[9, 11] (for the appropriate w), we observed several orders of magnitude difference for n > 1 whereas the dominating 
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FIG. 5: Value of $3 as a function of r, for uj = 0.86. The first panel shows the region near the origin, the second one the region 
where the oscillations begin to dominate and the third one the region around the outer matching radius i?max. The position of 
-Rmax is denoted by the circle. 
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FIG. 6: The first panel shows the total mass of the system, the second one gives the value of $i at the origin and the third 
one the transition radius defined as the smallest radius for which $3 = 0. All quantities are computed as a function of uj, in 
the high resolution case. 
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mode n = 1 is in relatively good agreement. We have currently no explanation for this discrepancy but feel rather 
confident with our results, given the various tests exhibited in Sec. IV B. As far as the metric fields are concerned, a 
direct comparison with [9, 11] is not easy because of the use of different systems of coordinates. Nevertheless, we are 
puzzled by the fact that the metric fields shown in [9, 11], at least for n > 0, do not seem to be even at the origin, as 
they should. It might be an effect of the resolution of the plots however. 
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FIG. 7: Values of the various modes for $ (first panel) and the metric fields A and B (second and third panels), for the 
oscillaton -with maximum mass, corresponding to a; = Wmin ~ 0.8608. 



The numerical values of various quantities are given in Table II. 
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TABLE II: Values of <I>i at the origin, the amplitude of the dominating tail C3 and the phase at which the mininmm is attained 
S3 . The values are given as a function of u and come from the computation with the highest resolution. 



OJ 


Mass 


$1 (r = 0) 


C3 


S-i 


0.83 


0.59446 


0.70326 


6.5894 ■ 10-'' 


+0.58 


0.84 


0.60036 


0.64172 


2.9848 


io-« 


+0.15 


0.85 


0.60399 


0.58468 


1.2793 


10"® 


-0.28 


0.855 


0.60495 


0.55760 


8.1692 


10-^ 


-0.50 


0.8575 


0.60522 


0.54438 


6.4825 


IQ-'^ 


-0.61 


0.85875 


0.60530 


0.53785 


5.7641 


10"^ 


-0.67 


0.86 


0.60535 


0.53137 


5.1183 


10-^ 


-0.72 


0.860625 


0.60535 


0.52815 


4.8210 


lo-'^ 


-0.75 


0.86125 


0.60535 


0.52494 


4.5392 


lO"'^ 


-0.78 


0.861875 


0.60534 


0.52174 


4.2724 


10-^ 


-0.81 


0.8625 


0.60532 


0.51855 


4.0201 


10"^ 


-0.83 


0.865 


0.60515 


0.50593 


3.1399 


lO"'^ 


-0.94 


0.87 


0.60437 


0.48122 


1.8815 


10-^ 


-1.17 


0.88 


0.60093 


0.43380 


6.2311 


10"* 


+1.53 


0.89 


0.59484 


0.38878 


1.8127 


io-» 


+1.07 


0.9 


0.58588 


0.34589 


4.4800 


lo-'^ 


+0.62 


0.91 


0.57371 


0.30493 


8.9842 • 10"^" 


+0.16 


0.92 


0.55792 


0.26570 


1.3691 • 10^"' 


-0.28 


0.93 


0.53796 


0.22805 


1.4368 • 10^" 


-0.72 


0.94 


0.51310 


0.19187 


8.8780 • 10"^^ 


-1.15 



V. SMALL-AMPLITUDE EXPANSION 



Review of the formalism 



Oscillatons in the small amplitude lim.it can be well described by an asymptotic expansion in terms of a small- 
amplitude parameter [13]. At leading order, configurations with any frequency can be characterized by the localized 
solution of a pair of ordinary differential equations. Since we wish to compare the numerical results in Sec. Ill with 
those obtained by the small-amplitude expansion we give a short review of the main definitions and results in [13]. As 
the amplitude of oscillatons decreases, the geometry approaches the Minkowski metric, and the frequency approaches 
the mass limit from below, which has been set to 1 in our case. At the same time, the size of oscillatons grows without 
limit, as is also indicated by the asymptotic behavior (16) of the leading mode of the scalar field. This motivates the 
introduction of the rescaled radial coordinate 

p = er . (26) 

The relation between the the small-amplitude parameter e and the fundamental frequency uj can be chosen as 

£2 = l-cj2. (27) 

The field $ and the metric functions are then expanded in even powers of e as 



J2k 



k=l 

00 

A=l + ^e2fe^2fc, 
fc=i 

00 



(28) 
(29) 
(30) 



fc=i 



where we used tilde notation to distinguish the coefficients of the e expansion from the Fourier components defined 
in (9)-(ll). To leading order the field $ turns out to be proportional to e^ . From the assumption that the functions 
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remain bounded as time passes follows that the configuration has to be periodic, and the time dependence of the 
coefficients can be integrated out. Writing out the first couple of orders, 



$ = e'^p2 cos(cjt) + £^p4 cos{ujt) + e^pQ cos{u}t) + e 




cos(3cjt) + C'(£®), 



A=l + e^a2 + e 



,(") 



,(2 



a^"' + a^"^ cos(2a;i) +0{e'') 



B = l- e^a2 + e'^ 
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vi 



cos(2a;i) 



0(e« 



(31) 
(32) 
(33) 



where p2, P4, Pe, 02, a 



(0) „(2) 



and 64 are functions of the radial coordinate p. These can be obtained order by order 



by solving ordinary differential equations arising in the small-amplitude expansion. The leading order configuration 
is determined by the time-independent Schrodinger-Newton equations 



(34) 
(35) 



d'^S 2 dS ^^ _ 
rfp2 p dp 


= 


d^s 2 ds ^2 
dp^ p dp 


= 



where the functions s and S are defined by 

s = -l-a2, S' = p2\/2. 
The total mass, M — 2ro, of the configuration can be expanded as 

M = eMi -I- e^M2 + O (e^) , 



(36) 



(37) 



where the numerical values of the constants are Mi = 1.75266 and M2 — —2.11742. 

The small-amplitude expansion gives exponentially localized configurations to all orders, and consequently it cannot 
describe the oscillating tail responsible for the mass loss of oscillatons. Although the small-amplitude expansion is an 
asymptotic one, it gives a useful representation of the core region. The dominant radiating mode (17) in $3 can be 
calculated to leading order by the extension of the Fourier mode equations into the complex r plane, and employing 
Borel summation in a region around the resulting pole [13], giving 



Cs 



kirQ 



■ exp 



\/8Q' 



(38) 



where Q is the distance of the pole from the real axis of the solution of the Schrodinger-Newton equations, and k is 
a constant. Numerically Q ~ 3.97736, k = 0.301 and 



C. 



3.761 / 11.2497\ 



(39) 



B. Comparison with the numerical solution of the Fourier mode equation 



In Fig. 8 we again plot the highest resolution tail amplitude data of $3 from Fig. 4, this time as a function of 



the small-amplitude parameter e = yl 
following empirical formula 



- UJ^ 



The theoretical curve from (39) goes well below the data points. The 



^emp. ^ 3J61 ^^ 



2\ 16.63 

£ ) exp 



11.2497 



(1-0.2990 £2 



(40) 



gives an excellent fit to the C3 data points, moreover, for small e it approaches the theoretical result (39). A natural 
interpretation of this result is that the number Q in the exponent of (38), describing the distance of the pole from 
the real axis, has a polynomial e dependence, and the constant k also changes for large e. 

The ratio of the numerically and theoretically obtained tail amplitudes is plotted in Fig. 9. The numerical value 
is several hundred times bigger for large e but the ratio decreases exponentially when e gets smaller. Even for the 



17 



0.0001 



1e-06 




1e-14 



1e-16 



Numerical o 

e-expansion 

« emp. 

'^3 



0.3 0.35 0.4 0.45 0.5 0.55 0.6 

e 



FIG. 8: Comparison of the numerical and theoretical values of C3. A very good empirical fit, given by (40) is also shown. 
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FIG. 9; Ratio of the numerical and theoretical values of C3, together with a fit of an exponential curve through the data points. 



smallest e where the tail could be calculated, with amplitude of about 10^^^, the difference is about twentyfold. 
Although because of numerical errors we cannot calculate the tail for smaller e, it appears convincing that the two 
method would converge for even smaller amplitudes. It can be seen in Fig. 9 that the fit of the exponential curve 



^£-exp. 



= 0.08045 cxp(15.870£) 



(41) 



to the ratio of the numerical and analytical result gives a very good approximation for 0.3 < e < 0.6, but is clearly 
not appropriate for smaller e values, since in the e — )■ limit it does not approach the value 1. 

The reason for the large difference between the theoretical and numerical tail amplitudes can be better understood 
by inspecting how precisely the core region is described by the small-amplitude expansion when e goes above 0.34. To 
do this, we compare the Fourier components of $ at the origin. From (31) it follows that $1 = £^P2 + £^P4 + 0{e^). 
Substituting the numerical values of p2 and p4, for the central value of the leading Fourier component we obtain 

^e-^oxp. ^ 2.44461 e^ + 1.49305 e"^ + 0{e^) . (42) 

In Fig. 10 we plot the relative difference of the numerical and the analytical results for $ic. Here we can include 
numerical values even for e < 0.34, since the core region can be reliably calculated by numerically solving the Fourier 
mode equations even if the tail amplitude in $3 goes below the numerical noise. We can see that the first two 
terms given by the small-amplitude expansion (42) give excellent approximation to the central amplitude for low and 
moderate e values. For the largest e the difference between the numerical and theoretical values grows to about 20 
percent. At £ = 0.34 it decreases to about 2 percent. 
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FIG. 10: Relative difference of the numerically calculated central value of the "l>i Fourier mode and the one calculated by the 
£ expansion method, using the e^ and e** order approximations from (42). To show the tendency of the points two curves 
proportional to e^ and e* are also plotted. 



The error of the s expansion for s > 0.34 is considerably larger for the next Fourier component $3, which component 
also gives the dominant contribution to the oscillating tail responsible for the energy loss. Using (31) the small- 
amplitude expansion gives the central value 



nr- 



0.0343467£'^ + C'(e^) . (43) 

In Fig. 11 the relative difference of the numerical value of $3c with respect to that given by (43) is shown. Here, for 




0.01 



FIG. 11: Relative difference of the numerically calculated central value of the $3 Fourier mode and the one calculated by the 
e expansion method. To show the tendency a curve proportional to e^ is also plotted. 

the largest parameter value, at e = 0.5578, the numerically calculated $3c is about five times the one given by the 
small-amplitude expansion. For e ~ 0.34 the numerical value is still about 50% bigger. Since the result (38) is based 
on the assumption that the core is described by the small-amplitude expansion, this big difference in the central value 
of $3 makes the large error in the tail amplitude for e > 0.34 comprehensible. 

In Fig. 12 we plot the total mass of the oscillaton as a function of the e parameter. The first two terms in the 
small-amplitude expansion, given by (37), give a reasonable approximation, but the position and the value of the 
maximum is not extremely precise. This not so surprising because the maximum is reached for s ~ 0.5 which is not 
particularly small. It is possible to make a very good empirical fit of the form 



M° 



eMi + e^M2 + e^Ms + e^ M^ + e^M^ 



(44) 
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FIG. 12: Numerically and analytically obtained masses of configurations with various e parameters. A quite precise empirical 
fit given by (44) is also shown. 



where M\ and M^ was kept at the value M\ — 1.75266 and M^ = —2.11742 given by the small-amplitude expansion, 
and the fitted constants are M3 = -0.24723, M4 = 1.1749 and M5 = -4.1308. The numerically obtained value 
for the place of the maximum is Emax — 0.509, corresponding to the frequency Wmin ~ 0.8608, and the maximal 
mass is M^ax = 0.60535. The importance of these considerations lies in the fact that oscillatons with e > emax, or 
equivalently, to < Wmin are unstable. The value given in [10, 11] for Wmin is 0.864 and for Mmax is 0.607. Reintroducing 
the scalar field mass, to, and expressing it in units currently used in Particle Physics, i.e. in eV/?nc^, we obtain the 
maximal oscillaton mass in kg-s: 



Mmax = 1.6085 X 102"kg 



eV 



mc^ 



(45) 



VI. MASS LOSS RATE 



In order to calculate the rate by which the mass M of the oscillaton decreases, we have to calculate the mass-energy 
carried out by the spherical scalar wave 



$ 



C3 



cos (Aar — 3ujt) , 



(46) 



where A3 — -\/9w^ — 1. Restoring the I/a/Stt factors from (3) into the scalar field, the mass-energy carried by the 
wave is 



dM 
~di 



d^ d$ 



2 dt dr 
Substituting (46), taking the large r limit, and averaging in time 

Using the empirical formula (40) for C3, for the mass- loss rate wc obtain 



dM 
IT 



-10.6l "^^t'' (l + ^')''"' 



exp 



22.4993 



(1- 0.2990 e^ 



(47) 



(48) 



(49) 



where lo = y/l — e^. Using the version of (48) with ui ~ 1, and substituting C3 from (39), we get the leading order 
small-amplitude result already obtained in [13], 



dM 30.0 / 22.4993\ 

^=-— ^"H — — 



(50) 
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Although for small e the empirical mass-loss rate formula (49) can be approximated by (50), for close to maximal 
amplitude oscillatons (49) gives significantly higher radiation loss. 

In order to be able to draw physical conclusions, we have to restore the scalar field mass m into the expressions. This 
can be done by substituting t — >■ mt, r — >■ mr into the equations (see e.g. [13]). Then wc have to substitute M — >■ mM, 
so the right hand sides of (37) and (44) will get a 1/m factor. The m factors drop out from (48), (49) and (50), so 
they remain valid for any scalar field mass m. The physical frequency of the oscillaton will he uj ~ muj ~ m^/l — s^. 

Taking (49) at Emax = 0.509 corresponding to the maximal mass M^ax = 0.60535 and dividing by M^ax we get the 
relative mass loss rate for the heaviest stable oscillaton: 



1 dM\ 



-5.917 X 10""r 



(51) 



where m is the scalar field mass in Planck units. This is about 14000 times larger than the theoretical estimation 
—4.3 X 10^^^ in [13] obtained from the leading order small-amplitude behavior. The large difference arises because 
the small-amplitude analysis underestimates the amplitude of the radiating tail by about a factor of 100 for e as large 
as £max = 0.509. 

It is natural to start with a maximal mass configuration with M = Mmax, and ask for the time period until the 
mass decreases to a certain part of its original value. Since the elapsed time t is inversely proportional to the scalar 
field mass m, in table III wc list the product tm. We give the elapsed time calculated first by using the small- 
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0.01 
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1400 
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10^1 
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^q30 
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10^" 
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1.99 • lO^-' 
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^q62 
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0.142 
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0.107 
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10*6 
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TABLE III: The time necessary for an initially maximal mass oscillaton to lose the given part of its mass. The value of tm is 
given in Planck units, and also when the time is measured in years and the scalar mass in eV/c? units. The values calculated 
by the small-amplitude expansion method and by the numerical solution of the Fourier mode equations are also given. 

amplitude expansion results and then by the more precise Fourier mode decomposition method. In the first case we 
calculate -^j- from (50), and approximate the e dependence of the mass by (37). In order to obtain the more precise 
result we calculate the mass- loss rate from the empirical formula (49), and approximate the mass function by (44). 
Unfortunately, a sign mistake was made in [13] when substituting the numerical value of M^^> during the calculation 
of the numbers in Table VII and VIII, which resulted in smaller radiation rate for large amplitude oscillatons than 
the correct rate. This is also the reason why we also present the numbers obtained by the e expansion here. 

Next we address the question that how much of its mass an initially maximal mass oscillaton loses during a time 
period corresponding to the age of the universe, which we take to be 1.37- 10^" years. In Table IV we list the resulting 
oscillaton masses in units of solar masses {Mq) as a function of the scalar field mass in eV/c^ units. 



VII. CONCLUSIONS 



We have presented a rather detailed numerical study of the structure of spherically symmetric, time periodic 
oscillaton solutions of the Einstein-Klein-Gordon equations. We solved the equations by using a two-dimensional 
spectral method for both the radial coordinate and time. The use of spectral methods enabled us to obtain very 
precise solutions at a moderate computational cost. 

For the first time we have succeeded in computing the amplitude of the standing wave tail of the time periodic 
oscillatons. The amplitude of those tails has been found to be very small indeed as compared to the central amplitude 
of an oscillaton. This implies that truly localized, time-periodic, asymptotically flat oscillatons do not exist, rather 
oscillatons of finite mass created by physical processes continuously lose some of their mass due to scalar radiation. It 
should be noted, however, that since the radiation rate of oscillatons decreases sufficiently rapidly as their total mass 
decreases, oscillatons cannot radiate away their mass in a finite time. 
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eV 


^-expansion 


Fourier 


expansion 


S^max £■ 


M 


M„ax 


^max ^ 


M 


-M„ax-M 
M„ax 


^Q-35 


2.31 • 10-1° 


8.20 


10^^ 


2.91 • lO-i'* 


2.37- 10"* 


8.09 


10^4 


3.75 • lO-i"' 


10-30 


7.31 ■ 10-« 


8.20 


10" 


2.91 ■ lO-i* 


7.50 ■ 10-" 


8.09 


10i« 


3.74 • 10-"' 


jQ-25 


2.31 • lO-'' 


8.20 


10" 


2.90 ■ 10-'-* 


2.18- 10-3 


8.09 


lOi* 


3.16 • 10-'' 


10-2° 


0.00621 


8.20 


10^ 


2.09 ■ 10-1 


0.0544 


7.94 


10« 


0.0183 


10-"' 


0.0883 


7.87 


10^ 


0.0400 


0.126 


7.36 


10* 


0.0896 


10-"' 


0.169 


7.06 


10-1 


0.139 


0.183 


6.65 


10-1 


0.178 


10-5 


0.226 


6.24 


10-" 


0.238 


0.227 


5.96 


10-" 


0.263 


1 


0.267 


5.56 


10-" 


0.322 


0.262 


5.36 


10-11 


0.337 


10^ 


0.298 


4.98 


10-1" 


0.392 


0.289 


4.85 


10-1" 


0.401 


10"^ 


0.323 


4.51 


10-21 


0.450 


0.311 


4.41 


10-21 


0.454 


10^5 


0.342 


4.11 


10-2" 


0.499 


0.329 


4.04 


10-2" 


0.500 



TABLE IV: Mass M of an initially maximal mass oscillaton after a period corresponding to the age of the universe for various 
scalar field masses. The decrease in e from £max, and the relative mass change (Mmax — M)/Mmax is also given. At the 
small-amplitude expansion Emax = 0.525 is used, while the mode decomposition value is £max = 0.509. 

Using the precise numerical results we have derived a semi-empirical mass loss formula of "physical" oscillatons 
of finite mass. The results show that the previous computations of the mass loss rates in the small amplitude 
limit underestimated the true rate by several orders of magnitude for larger amplitude oscillatons. Nevertheless 
the qualitative picture of an oscillon as a lump losing its mass extremely slowly prevails making these objects of 
physical interest, such as dark matter candidates. The agreement with analytical results obtained in the limit of 
small amplitudes is very satisfactory as far as the structure of the core is concerned. Concerning the amplitude of 
the oscillatory tail it is more difficult to do a precise comparison. Indeed, the tail is very small and can be accurately 
extracted from the numerical simulations only when the amplitude of the oscillatons is not small, so in the region 
where the analytical approximation is expected to fail. The numerical and analytical results are coherent with each 
other. 

We have also computed the value of the maximal mass which an oscillaton may have, together with the corresponding 
value of frequency. 
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